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Abstract 


This research is involved with the implementations of advanced computational schemes 
based on large eddy simulations (LES) and direct numerical simulations (DNS) to study the 
phenomenon of mixing and its coupling with chemical reactions in compressible turbulent 
flows. In the efforts related to LES, we have initiated a research program to extend the 
present capabilities of this methods for the treatment of chemically reacting flows, whereas 
in the DNS efforts, we have focused on detail investigations of the effects of compressibility, 
heat release and non-equilibrium kinetics modelings in high speed reacting flows. Our 
efforts to date, have been primarily focused on the simulations of simple flows, namely 
homogeneous compressible flows and temporally developing high speed mixing layers. 

This report provides a summary of our accomplishments during the first six months of the 
research supported under Grant NAG 11122 sponsored by NASA Langley Research Center 
administrated by Dr. J. Phil Drummond of the Theoretical Flow Physics Branch. 
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Summary of Accomplishments to Date 


A review of the recent computational work on internal supersonic reacting flows indicates 
the need for further development in both the methodology and implementation of DNS 
and LES for the treatment of the turbulent reacting flows. The results of previous 
and ongoing works at NASA Langley toward the developments of advanced numerical 
algorithms for the simulations of such flows have been very successful, and we have used 
these computational tools as a means of understanding some of the interesting intricate 
flow dynamics in conjunction with the implementations of DNS and LES in supersonic 
combustion. 

In this work, we have initiated a rather broad research program in development and im- 
plementations of both DNS and LES for the treatments of compressible reacting flows. 
Below, a summary of our efforts during the first 6 months of this research is provided and 
detail descriptions are furnished in Appendix I and Appendix II. 

Large Eddy Simulations of Compressible Reacting Flows: 

Our major goal in this effort is to initiate a program to extend the present capabilities 
of LES for the treatment of chemically reacting flows. In the efforts to date, we have 
been primarily concerned with a 'priori analysis of subgrid fluctuation in a compressible 
homogeneous flows. This analysis is mainly involved with constructing the shapes of the 
simulated PDF’s within the subgrid. This work is currently in progress, and Appendix I 
provides a summary of our state of progress. 

Direct Numerical Simulations of High Speed Reacting Mixing Lavers: 

This work is involved with the assessment of the roles of compressibility and heat release, 
and the non-equilibrium effects of chemical kinetics on the compostional structure of the 
flame in high speed reacting mixing layers. This work has been mainly motivated due 
to the iterest of NASA in understanding the global and detail mechanisms of mixing and 
chemical reactions in high speed flows. Appendix II provides a detail description of our 
achievements to date. 
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On DNS and LES of Reacting Compressible Homogeneous Turbulence 


Cyrus K. Madnia and Peyman Givi 
Department of Mechanical and Aerospace Engineering 

SUNY at Buffalo 
Buffalo, NY 14260 

ABSTRACT 

Direct numerical simulations (DNS) are performed for two and three dimen- 
sional compressible reacting turbulent flows. A range of initial fluc- 

tuating Mach number was considered to study the effect of compressibility 
on mixing and chemical reaction. A passive second order chemical reac- 
tion of a type A + B — ► Products was considered. As am initial effort, 
it is assumed that the chemistry is infinitely fast (i.e. Damkohler Num- 
ber — > oo) ; therefore a flame sheet approximation is employed. With this 
approximation, the transport of an inert scalar quantity is sufficient 
to portray the statistical behavior of the species field. A spectral- 
collocation algorithm based on Fourier expansion function 1 is employed in 
the numerical simulations. The fluid obeys the Navier-Stokes equations 
for an ideal gas. Periodic boundary conditions are employed in all di- 
rections for both two dimensional and three dimensional simulations. A 
256 x 256 uniform collocation grids was used for two dimensional calcula- 
tions. The number of collocation points for the three dimensional simula- 
tions was 128 x 128 x 128. 

For the purpose of flow visualization, two dimensional contour plots of 
density are presented in Figs. 1 and 2. Figure 1 represents a typical 

density contour plot for a low compressible case and Fig. 2 corresponds 
to a high compressible case. An interesting feature displayed in Fig. 

2 is the presence of the shocklets. A two dimensional contour plot of 
species concentration for a high compressible case is shown in Fig. 3. 

The normalized time corresponding to this case is such that the effect of 
initial conditions is no longer substantial. 

The DNS data base obtained is used to study the behavior of the probabil- 
ity density functions (pdf’s) of scalar properties within the subgrid. 

The use of pdf methods in stochastic description of reacting flows has 
proven valuable in Reynolds averaging turbulence modeling 2,3 and it is ex- 
pected that their implementation for subgrid closures in LES of reacting 
flows would also be beneficial. After the generation of the data base on 
the fine grid, the results are statistically analyzed within an ensemble 
of these grids to describe the large scale conduct on the coarse grid 3 ' 4,5 . 
The ratio of the mesh spacings (resolution) provided by the coarse and the 
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fine grids is a measure of the size of the filter which would be used in 
LES. This analysis is done at intermediate computational times, at which 
the influence of the initial conditions (of the chemical field) is not 
substantial. The results of these a priori evaluations indicate that the 
pdf of the inert scalar within the subgrid (of various filter widths) re- 
sembles that of a Gaussian distribution. This had been already surmised 
in incompressible flow simulations, as previous DNS results 6,7 had sug- 
ij 3 t'd. in present simulations this behavior is observed both 

in incompressible flows and in compressible flows dominated with shock- 
lets. This observation is somewhat useful since it suggests that in 
subgrid modeling of an inert scalar property, the information on the first 
two moments of the variable is enough to parametrize the pdf within the 
subgrid. 

With the knowledge gained to date, it is anticipated that the approach 
based on pdf parametrization based on its first two moments may prove ser- 
viceable for turbulent combustion simulations. The approach based on 
the solution of a transport equation for pdf, however, may not be prac- 
tical at this stage. An estimate of computational requirements indi- 
cates that the cost associated with LES (a semi-deterministic solution of 
large scale with a probabilistic description of small scale by solving a 
pdf transport equation) is of the same order as that of DNS on the fine 
grids, unless the ratio of the fine to coarse grid is large. 

Our ongoing investigation is concerned with investigating the effects of 
finite Damkohler number, which is most appropriate for pdf modeling, and 
also on including the influence of the heat release. The statistical 
analyses are also being done for different flow types and for various fil- 
ter widths. 
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Figure 1. Contour plot of density for low compressible case. 
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Figure 2. Contour plot of density for high compressible case. 
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Figure 3. 


Contour plot of species A concentration for high compressible case. 
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Abstract 


The effects of compressibility, chemical reaction exothermicity and non-equilibrium 
chemical modeling in a reacting plane mixing layer were investigated by means of two 
dimensional direct numerical simulations. The chemical reaction was irreveisible and 
second order of the type A + B -*• Products + Heat. The general governing fluid 
equations of a compressible reacting flow field were solved by means of high order finite 
difference methods. Physical effects were then determined by examining the response' 
of the mixing layer to variation of the relevant non-dimensionalized parameters. 

The simulations show that increased compressibility generally results in a sup- 
pressed mixing, and consequently a reduced chemical reaction conversion rate. Re- 
action heat release was found to enhance mixing at the initial stages of the layer's 
growth, but had a stabilizing effect at later times. The increased stability manifested 
itself in the suppression or delay of the formation of large coherent structures within 
the flow. 

Calculations were performed for a constant rate chemical kinetics model and an 
Arrhenius type kinetics prototype. The choice of the model was shown to have an 
effect on the development of the flow. The Arrhenius model caused a greater temper- 
ature increase due to reaction than the constant kinetics model. This had the same 
effect as increasing the exothermicity of the reaction. Localized flame quenching was 
also observed when the Zeldovich number was relatively large. 
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Chapter 1 


Introduction 


1.1 Motivation 

Research on the subject of compressible reacting mixing layers has been of high prior- 
ity in recent years [1], Much of this effort has been devoted to the development of high 
speed air breathing flight vehicles. This type of vehicle would, according to current 
proposals, use a supersonic combustion ramjet (scramjet) engine. In such an engine, 
fuel is injected into a high speed airflow. The mechanisms of mixing and combustion 
of this non-premixed, high speed, compressible flow is of great complexity. In a sim- 
plistic approach, the problem may be divided into two parts: the effect hydrodynamic 
phenomena have on combustion, and the effect the combustion processes have on the 
hydrodynamics. Even divided as stated, the problem is still notoriously difficult. As 
such, only a few aspects of such phenomena and effects can be investigated in a single 
study. 

In this investigation, it was intended to study some of the effects of coupling 
between mixing and reaction in compressible combusting systems. The system con- 
sidered was a two-dimensional planar mixing layer, and a computational approach 
was chosen. The computational tool used was direct numerical simulation based on 
higher order finite difference algorithms. This method was chosen because of its free- 
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dom from turbulence modeling requirements (known also as a “model-free simulation” 
[2]), and its ability to compute the details of non-linear physical phenomena. The 
computational approach is suitable for use as a tool for basic research or as a compli- 
ment to experimental studies. The primary advantage of direct numerical simulation 
is the capability of evaluating all pertinent statistics of the flow without resorting to 
turbulence closure models. Simulations such as this usually require computational 
resources that were not readily available previously to researchers. 


1.2 Previous Research 

Given the extent of research on high speed mixing layers recently, it is not possible 
to include a comprehensive review of the accomplishments in this research area in 
this work. Thorough surveys of the state of the art have been performed by Givi 
and Riley [1], Givi [2], and Drummond [3], and will not be repeated here. Instead, 
a summary of some previous work will be presented, with priority given to those 
with direct relevance to the present effort. In particular, the scope of this review 
is limited to describing some of the recent accomplishments in investigations of the 
effects of compressibility on turbulence and large scale structures in parallel shear 
flows, and the influence of heat release and nonequilibrium chemical reactions on the 
development of these structures. 

Brown and Roshko [4] found that the turbulent mixing layer is dominated by large 
scale coherent structures, or vortices. These structures convect at a nearly constant 
speed and tend to coalesce with neighboring vortices. The authors demonstrated that 
the reduction in mixing layer growth rate that had been observed in their experiments 
was due to the influence of compressibility effects, not density variations as had been 
thought previously. Ho and Huang [5] showed how the growth rate ot the mixing layer 
could be manipulated by perturbing the flow at a subharmonic of the most amplified 
frequency. This technique stimulates the merging of the vortices, thus accelerating the 
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growth of the layer. A thorough review of the effects of harmonic forcing techniques 
is available in [6]. 

Papamoschou and Roshko [7] continued these experiments to examine the efFect 
of compressibility on the spreading rate of a supersonic mixing layer. They found 
that it is useful to study the flow in a reference frame that travels with the flow at 
the same speed as an average large scale structure. .. parameter which •> •' n .ti p es 
the compressibility in the flow was proposed as the convective Mach number, M c \ 
defined as the Mach number of the flow with respect to the above mentioned frame 
of reference. A direct correlation was found between M c and the stability of the flow. 

Ragab and Wu [8] substantiated the use of the convective Mach number as a rel- 
evant compressibility parameter by analyzing linear instability waves in supersonic 
shear layers. They also determined the stabilizing efFect of the velocity and the tem- 
perature ratios between the two streams of the flow on the stability of the layer. It 
was found that there is a complex, non-linear relationship between the growth rate 
of the waves and the velocity ratio. It was also shown that at low Mach numbers, a 
temperature increase has a stabilizing effect, whereas at high Mach numbers the effect 
is to destabilize the flow. Lele [9] verified the results of Papamoschou and Roshko, 
by means of direct numerical simulation of a two-dimensional layer. He proposed an 
explanation of the compressibility stabilization efFect based on the inviscid vorticity 
equation. Also, the development of eddy shocklets in the flow was noted for M c > 0.7. 
These shocklets are formed as a result of locally supersonic regions that appear dur- 
ing vortex roll up or pairing, and remain attached to the vortices as the structures 
travel downstream. The effect of compressibility was further studied by Sandham and 
Reynolds [10] for both two and three dimensional mixing layers. The mixing layer 
growth rate was found to be reasonably predicted by linear stability analysis. Shock- 
lets were captured in two-dimensional simulations when M c > 0.7. It was also found 
that three-dimensional effects become significant at M c > 0.6, and become dominant 
at M c > 1 . However, no eddy shocklets were observed in three-dimensional simula- 
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tions. Elliott and Samimy [11] performed experiments to investigate high Reynolds 
number compressible flows. They found a reduction in the level of turbulence fluctu- 
ations as the convective Mach number is increased. 

The effects of an exothermic chemical reaction on fluid dynamics is another area of 
interest. Among recent computational efforts, McMurtry, et al. [12] performed direct 
numerical simulations of turbulent reacting mixing layers. They used an approximate 
set of equations that are valid for low Mach number flows. It was found that the heat 
liberated from a chemical reaction causes the layer to grow at a slower rate than that 
of a non-heat releasing flow. It was also shown that the magnitude of product formed 
and the amount of mass entrained into the vortical structures decrease as the intensity 
of heat release increases. These results agree with those obtained experimentally (e. 

g- [13])- 

Jackson and Grosch [14] performed a linear stability analysis on supersonic react- 
ing mixing layers. They found the existence of fast and slow stability modes, with the 
slow mode appearing only for flows with heat release. An increase in the amount of 
heat release was found to result in a reduction in the growth of the fast waves along 
with an increase in the growth of the slow waves. As the heat release is increased 
to large values, the slow mode becomes the most unstable. Thus, it was determined 
that the overall effect of increasing the heat release is to first stabilize the flow, then 
to destabilize it. 

A direct numerical simulation of a supersonic reacting mixing layer was performed 
by Sekar, et al. [15]. They found reductions in the convective speed, the growth rate, 
and the entrainment of the free stream flows with increase in the magnitude of the 
heat release. They also found that the reduction of turbulence fluctuations with heat 
release occurs in supersonic flows to a lesser extent than that in incompressible flows. 
Their final conclusion was that heat release may not have a significant influence on the 
structure of the flow. Therefore, they suggest that for investigations concerned with 
mixing effects it might be useful to concentrate on phenomena related to gas-dynamic 
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effects rather than on exothermicity. 

The phenomenon of flame extinction in non-premixed flames was the last area 
to be investigated. Although this has been the topic of theoretical and experimental 
study, direct numerical simulations of such flows have been somewhat limited. Recent 
reviews of some of the prevalent theories regarding the structure of turbulent non- 
premixed flames, as well as some of the experiments and numerical work in thct fi ;ld 
are provided by Bilger [16] and Peters [17]. A point that is made in these reviews 
is that although turbulent combustion modeling is very useful, there exists great 
uncertainty in the formulation of these models and in their use. This uncertainty is 
avoided by using “model free” simulations. Such simulations of non-premixed flames 
have been performed by Givi, et al. [18] and Givi and Jou [19]. However, these 
studies made certain limiting assumptions (i. e. constant density) that limit their 
applicability to low speed flows. 


1.3 Scope of Present Research 

The objective of this work is to examine the effects of compressibility and chemical re- 
action exothermicity on a reacting plane mixing layer. An examination is also made 
of the non-equilibrium effects of the chemical kinetics on the structure of a flame. 
These are accomplished by direct numerical simulation of an unsteady two dimen- 
sional layer. The governing equations are integrated via high order finite difference 
methods. Physical modeling is kept as simple as possible so that the physical effects 
described in the previous section can be isolated. This has an added advantage of 
saving considerable amounts of computational resources. 

The mixing layer is assumed to be “temporally developing” with periodic bound- 
ary conditions. The fluid is assumed to be calorically perfect and to have constant 
and identical thermodynamic parameters. A simplified one step, second order irre- 
versible reaction is used to describe the reactant conversion. Both constant rate and 
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Arrhenius type kinetics models are used. In most cases, the layer is perturbed only 
by numerical truncation and round off errors. For those cases where these are not 
sufficient to destabilize the flow, explicit harmonic forcing is added. 

Physical effects are studied by changing the appropriate nondimensionalized pa- 
rameters. Compressibility is represented by the convective Mach number, and re- 
action exothermicity is measured by the nondimensional value of the enthalpy of 
reaction. The chemical reaction is controlled by the magnitude of the Damkohler and 
the Zeldovich numbers. 

The problem to be solved is formulated in Chapter 2, where the governing equa- 
tions are presented and discretized into a vector form. The physical models used in the 
simulations are then presented, followed by a description of the numerical algorithms, 
and boundary conditions. Results from the simulations are presented in Chapter 
3, where the effects of compressibility, reaction exothermicity and non-equilibrium 
chemistry are discussed. Finally, a summary, conclusions, and recommendations for 
future work are presented in Chapter 4. 
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Chapter 2 


Numerical Formulation 


2.1 Governing Equations 

A two-dimensional compressible, reacting flow is governed by the continuity, momen- 
tum, energy, and species conservation equations coupled together with an equation 
of state. These are expressed as [20]: 

Continuity 



| + v.(,n=o 

(2.1) 

Momentum 

+ V . (pVV) = v-T + p£f,b, 

(2.2) 

Energy 

d( P E) 

dt 

+ V . (pVE) = v . (r . V) - V . q + p^fMV + Vi) 

i=i 

(2.3) 

Species continuity 

+ V • (pVfi) = Wi-V ■ (pfiVi) 

(2.4) 

Equation of state 

f 

p = pRT ^W 
1=1 iw * 

(2.5) 
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where 


_ (dui duj\ du k 

r_r 0 = %- + /‘( dl . + aii 

(2.6) 

and 

q=-kVT + pf^h t f t V t 
» = 1 

(2.7) 

The total energy is given by 


N * v 2 u 2 

EW-f + Ef 

i=l r ,=i z 

(2.8) 

and the enthalpy of species i is defined as 


hi = h°+ f T c Pt dT i = 1,2 ,...,iV 5 

J T r 

(2.9) 


2.2 Physical Modeling 

This section discusses the various models needed to describe molecular diffusion and 
chemical reaction. For both processes it is possible to give complicated and computa- 
tionally expensive models. However, to keep the computational cost at an affordable 
level, some simplifying assumptions are made. Since the goal of this research is to 
investigate the physics of reacting plane mixing layers in a general sense; models that 
are limited to specific reactions or particular species are not considered 

It is assumed that the diffusivities of each chemical species are the same. There- 
fore, the diffusion of a species into another is proportional to respective concentration 
gradients. The effect of this is to decouple the equations for the diffusion velocities, 
and results in a form of Fick’s law: 


Vij = - 


DdU_ 

/. dx } 


( 2 . 10 ) 


where is the diffusion velocity vector of the ith species in the j th coordinate 
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direction and D is the binary diffusion constant. The value of this constant is deter- 
mined by choosing an appropriate value of the Schmidt number, Sc, since D = 
Analogously, the mixture thermal conductivity is expressed by 


k = 


SeU 

Pr 


( 2 . 11 ) 


where Pr is the Prandtl number. 

In earlier versions of the computational methodology used in this study [21] a dif- 
fusion model based on kinetic theory was used. In that model the diffusion velocities 
are described by the solution of an equation of the type: 

Ns Y Y \7n 

vx, = E 1 i r L (v i -Vi) + (fi-x i )-^ + 

i=i D a p 

^XiXjfDrj D t ,\ VT 

U p ° y V Si Si ) t 

The solution ofEq. (2.12) requires solving a system of N s simultaneous equations, with 
N s representing the total number of species involved in the reaction. This solution 
is computationally intensive and requires a coupled system of equations at each grid 
point throughout the computational domain. This process can require as much CPU 
time as solving the Navier-Stokes equations for the convective velocities [22]. 

In this work all of the species within the flow were given identical thermodynamic 
properties. Also, the assumption of a calorically perfect gas was made. All external 
body forces 6, were assumed to be negligible. 

The chemical reaction in the flow is assumed to be of a simple, irreversible, second 
order type of the form 

A + B — ♦ Product + Heat (2.13) 

The reaction is characterized by the kinetics mechanism, which is given by the single 
step model of 

w = KjC a C b (2.14) 


p )=1 
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where Ca and Cb represent the concentrations of the reacting species and are assumed 
equal at the free streams, i.e. Ca «, = Cb «, = Coo- Kj is the reaction rate constant, 
and can be normalized to form the definition of the Damkohler number, Da: 


= KjC oo 

Doo / 6 w \q 


(2.15) 


In the present study two types of chemistry models were used; constant rate kinet- 
ics (i.e. constant K /) and an Arrhenius type model in which Kj varies with the 
temperature. This is written as 


K f = A } e TTfc 


(2.16) 


where A/ is the pre-exponential factor and Ze is the Zeldovich number, defined as 


Ze = 


E 

RToo 


(2.17) 


Here E is the activation energy and R is the universal gas constant. When the 
Arrhenius kinetic model is used, the pre-exponential factor Aj replaces Kj in the 
definition of Da (Eq. 2.15). 

Combustion exothermicity is measured by the energy liberated by the chemical 
reaction, AH 0 . The magnitude of this energy is parameterized by a non-dimensional 
heat release parameter Ce, defined by: 


Ce = 


-A H° 


CvToo 


(2.18) 


Thus, Ce = 0 corresponds to a non-heat releasing chemical reaction. 
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2.3 Numerical Solution of the Governing Equa- 
tions 

The next step in the formulation is the discretization and the integration of Eqs. 
(2-1)— (2.4). In a vector form these equations are expressed by: 


dV dF (U) dG(U) __ 

dt dx + dy 


where U is the dependent variable vector, 


(2.19) 


U = < 


P 

pu 

pv 

P E 

Pfi 


( 2 . 20 ) 


F and G contain the diffusive and convective flux vectors in the x and y directions 
respectively, 

( 

pu 

puu — (J x 

puv Ty X | ( 2 . 21 ) 

( pE - a x )u- T xy v + q x 
puf, + puj, 


F = l 


and 


G= l 


pv 

pUV - T X y 

PW — (Jy 

(pE Ty X U -|- Qy 

PV ft + pVifi 


( 2 . 22 ) 
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a x = t ix and a y = T yy . Finally, H is the source vector: 

0 

pEifAx 

P }Ay 

pHi / u .(v > V) 

Wi 

For the purpose of numerical discretization, it is convenient to map these equations 
from the physical domain into an appropriate computational space. In the simu- 
lations performed here the grids are highly compressed in the transverse direction 
of the flow, with maximum compression along the region of maximum shear. This 
compression provides a sufficiently fine resolution in the area of large velocity and 
concentration gradients. A detailed explanation of the grid generation routine and 
the transformation process may be found in [21]. 

Two numerical schemes were utilized to integrate the governing equations; an 
algorithm proposed by Gottlieb and Turkel [23] and a compact parameter scheme 
developed by Carpenter [24]. Both methods are second order accurate in time and 
fourth order accurate in space. The two algorithms are dissipative, allowing a more 
accurate treatment of sharp gradients compared with non-dissipative methods. The 
main advantage of these methods is their capability in capturing shocks and reaction 
zones. 

The Gottlieb-Turkel scheme is a variant of the well known MacCormack predictor- 
corrector method [25], For Eq. (2.19), it is implemented as: 

Predictor: 
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Corrector: 


U ■ ■ = U*, - f 


7F*,j - 8F*_j ; + F*_ 2 j] - f [7G-J - 8G*j_ 1 + GV_ a ] + AtH*, 


U ^ 1 = - fu n + U**l 

hj O [ ' MJ 


(2.25) 

(2.26) 


where \ x = ^ and \ y = The CFL (A r or A y ) condition for stability requires 
CFL < |. A disadvantage of the method is the need for a use of a five point stencil, 
thus precluding its use on the gridpoints next to the boundaries of the computational 
domain. For this reason, and in an effort to improve the accuracy, a family of dissi- 
pative compact parameter schemes (DCPS) are also considered. For Eq. ( 2.19), the 
DCPS takes the form of: 

Predictor: 


DSU tj + S I + D S? +1J = 


„ _ (B - X)F?. U + CF*, + (fl + 4)F? tu 


Ax 


PT „ i t» i nr - {B ~ a)g 6 >- + CG "> + 

+ ± i,i + 1 - ^ 


= U", - A((S", + T", - H",) 


Corrector: 


OS’-U + s h + DS Ui, = 


. ~(B + A)Fr_ u - CF-, - (B - A) 


Ax 


nr T . | n T . -(B + A)Qv_, - CG', - (B - A)Gv„ 

"6 A .,j + A ' A .,J+1 - 

U-^UV-AUSV+TV-HV) 

u ",r = 5 l u ”> + u ;.;l 


(2.27) 

(2.28) 

(2.29) 

(2.30) 

(2.31) 

(2.32) 

(2.33) 


In Eqs. (2.27)-(2.32) A = C = —25, and D = A [24]. S"^ and S*^ are the numer- 
ical values of the derivative of F, and T” ; - and T*j are the numerical values of the 
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derivative of G. S and T must be calculated implicitly at each of the predictor and 
corrector stages by inverting a tri-diagonal matrix. If a flow with periodic boundary 
conditions are considered, a periodic tri-diagonal matrix must be inverted. This repre- 
sents a one parameter family of methods in B. The value of B at which the maximum 
allowable CFL occurs is B = The CFL condition in this case is To preserve 
(At 2 , Ax 4 , Ay 4 ) accuracy, the predictor-corrector sequence musi oe switched at each 
time step (i.e. Forward/Forward-Backward/Backward, then Backward/ Backward- 
Forward/ For ward). This cycling procedure has an added advantage of dampening 
an instability that may occur when either scheme is formulated in two (or more) 
dimensions. 

The advantages of these higher order algorithms in comparison with some of the 
conventionally utilized lower order difference methods were demonstrated by com- 
paring results generated by all these schemes to the results of some test problems 
with known analytical solutions. The lower order discretization schemes are based on 
a first order upwinding method [26] and the second order MacCormack scheme [25] 
In this demonstration, a linear wave equation was considered, and comparisons were 
made for both one and two-dimensional cases. In the former, the linear advection of 
a square wave concentration distribution was considered, and in the latter the solid 
body rotation of a sharp gradient scalar field was investigated. With the absence of 
diffusion, the scalar field retains its initial shape in both cases; providing an effective 
means of evaluating the discretization routines. 

The results indicated that the Gottlieb-Turkel method and the DCPS provide a 
substantial improvement over the first order upwind and second order MacCormack 
schemes in that the magnitudes of both the truncation and the phase errors are sub- 
stantially reduced. Also, the DCPS method resulted in slightly lower phase error than 
generated by the Gottlieb-Turkel scheme. However, the computational requirements 
for the DCPS was about 25% more than that required by the Gottlieb-Turkel algo- 
rithm. Similar tests were previously performed by Carpenter [24], who utilized the 
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MacCormack, Gottlieb-Turkel, and DCPS methods to calculate the growth rates and 
the characteristic frequencies for a temporally developing compressible mixing layer, 
and compared the results to those given by a spectral linear stability method. In this 
comparison, the DCPS was found to be twice as accurate as the MacCormack and 
the Gottlieb-Turkel methods. 

In view of these comparisons, it was decided to select the Gottlieb- Tuikci and 
DCPS algorithms in favor of other alternatives. This decision was made mainly to 
keep the numerical truncation errors at most of order 0(At 2 , Ax 4 , Ay 4 ). In the 
subsequent chapters, both of these methods are used interchangeably with consider- 
ation to available computational resources. However, only one method was used to 
describe each physical phenomenon. Namely, the Gottlieb-Turkel method was used 
in the calculations discussed in Sections 3.1 and 3.2, whereas the DCPS was used in 
those presented in Section 3.3. 

2.4 Initial and Boundary Conditions 

One of the primary assumptions made in these simulations is that the mixing layer is 
temporally developing. That is, -the reference frame of the simulations is defined to 
be traveling along with the average velocity of the flow. The advantages of this ap- 
proximation are twofold. First, with the temporal assumption the inflow and outflow 
boundary conditions can be assumed periodic. This removes the difficult problem of 
specifying the boundary conditions. Second, the temporal assumption means that 
only a relatively small region of the flow is being simulated. This region is then 
followed in a Lagrangian sense as time progresses. This results in considerable com- 
putational savings, namely in CPU time and memory allocation. These savings can 
then be used to simulate the flow in a greater detail. The primary disadvantage of a 
temporal approximation is that asymmetric effects in the flow can not be captured 
[27]. These effects are not significant in the scope of the present research. 
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A representation of the initial flow field is shown in Fig. 1. The flow on the top 
stream moves towards the right with a free-stream streamwise velocity of Uoo > and 
the bottom stream moves to the left with the same speed ( — f/oo). The fiowfield 
is initialized with a hyperbolic tangent streamwise velocity profile with a specified 
initial vorticity thickness (<5 w |o). There is no initial fluid motion in the transverse 
(K) direction. The pressure was initially assumed to be constant throughout the flow 
field. Depending on the problem simulated, the temperature was either assumed to 
be initially constant or had an initial distribution of the form: 

T = 7^(1 +4e- 1000v ‘) (2.34) 

where y * is the normalized spatial coordinate in the transverse direction, defined as 
y* = (y — 0.5y ma x)/ymax- Reactant A covers the top half of the physical domain, 
and reactant B covers the bottom half. It was assumed that the upper and lower 
walls were far enough away from the mixing layer that free stream conditions could 
be imposed. 

In most cases, no explicit forcing was added to the base flow. The simulations 
relied on numerical truncation errors to provide perturbations to the layer to trigger 
formation of coherent vortices. For most cases, this was sufficient. However, in some 
cases where the physical effects had a stabilizing effect on the flow, harmonic forcing 
was explicitly added. The forcing used was that determined by a linear stability 
analysis of a temporally developing incompressible layer with the same initial velocity 
distribution as that employed here [28]. No attempts were made to find the most 
unstable modes for the compressible flow. This is justified in view of the fact that 
the study is focused on investigating the effects of large scale structures once they are 
formed, not how these structures are most rapidly generated. 
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Chapter 3 


Presentation of Results 


The results of the numerical simulations will be presented in three forms: integral 
or global representations, statistical sampling, and flow visualization. A global rep- 
resentation of how the state of the flow changes with time is given by considering 
the temporal variation of integral parameters, for example the vorticity thickness or 
the total amount of product formed. The second form of presentation of result is 
performed by examining the cross-stream variations of the statistics of the relevant 
variables. With the approximation of temporal evolution, the flow is considered ho- 
mogeneous in the streamwise direction X , and the statistical information is obtained 
by sampling the data in this direction. In this way, the ensemble mean and mean 
square of an arbitrary transport variable <p are obtained by 

(*(*')>-^5>(n (3-D 

iV 1 = 1 

<*>°<n) = ^£te -<¥>>)’ (3.2) 

^ 1=1 

where ( ) denotes an ensemble average, the subscript i indicates the grid index, and 
N is the total number of grid points in the X direction. 

The final method of result presentation, flow visualization, is provided by present- 
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ing the data in the form of contour plots. This type of presentation is an effective 
method of visualizing the flow. A typical contour plot of the vorticity is shown in 
Fig. 2. This figure shows that the computational grid is compressed in the transverse 
direction, with maximum compression occurring along the centerline of the vortex. 
Only a magnified portion of the entire computational grid is presented in this figure 
to highlight the details. All of the contour plots that follow will be presented on a 
uniform, or physical grid. The fluctuations near the outer portions of the figure may 
be attributed to numerical noise. Although the amplitude of the fluctuations is small 
compared to the physical values, they are still displayed due to deficiencies in the 
graphics software used. 

All parameters are normalized when appropriate by initial or free stream condi- 
tions. Time is normalized by 


t 

k/UZ 


(3.3) 


where Iq is the physical size of the computational box. 


3.1 Compressibility Effects 

The effect of compressibility on a mixing layer was studied by varying the convective 
Mach number, M c , while keeping all other parameters constant in a non-reacting 
layer. The flow was examined for M c = 0.2, 0.4, 0.8, and 1.2. The grid resolution 
of the M c = 0.2 case consisted of 127 x 127 grid points, and for the other cases it 
consisted of 256 x 256 grid points. The increased resolution was necessary to resolve 
the strong gradients that exist at high compressibility. No explicit forcing was added 
to the flow. The vorticity thickness vs. normalized time for different values of M c is 
given in Fig. 3. This figure shows that the growth rate of the vorticity thickness is 
clearly decreased with increased compressibility. Compressibility also affects the time 
needed for the layer to roll up into a vortex. The onset of roll-up is signified by a 
large jump in the vorticity thickness. The compressibility also affects the size of the 
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vortex formed, with the thickness of M c = 0.2 layer being about 50% larger than the 
M c = 0.8 vortex. 

The profiles of the normalized average streamwise component of the velocity is 
shown in Figs. 4-6. The most significant feature portrayed by these figures is the 
steepness of the mean velocity profiles at high convective Mach numbers. This shows 
a sharper velocity gradient across the layer, implying a lesser rate of mixing. As time 
progresses, the slopes of the profiles begin to decrease. The suppression of turbulence 
fluctuations with compressibility is shown by the profiles of the mean square of the 
fluctuating velocity. The transverse variation of the fluctuation is shown in Fig. 7 for 
time t* = 6 and in Fig. 8 for time f* = 8. A marked decrease in the amplitude of 
the fluctuations can easily be observed. For t * = 6, the fluctuations for the higher 
compressibility cases are almost negligible compared to those at low compressibility. 
At t* = 8, the fluctuations for M c = 0.8 and M c = 1.2 are evident, although at lower 
amplitudes than those in lower compressibility cases. The reduction of turbulence 
fluctuations is further illustrated in Figs. 9 and 10, which show Reynolds stress profiles 
in the mixing layer at times t* = 6 and t * = 8. 

The pressure response for the layer to the formation of large scale structures is 
shown in Figs. 11 through 16 for M c = 0.2 through M c — 1.2. The regions of pressure 
minima occur at the location of vortex cores. Similarly, pressure minima are located 
at the braids between the large scale structures, at the point of minimum vorticity. 
This is evident by a comparison between Figs. 11 and 25. Another interesting 
feature of the increased compressibility is shown by examining the plots of pressure 
contours at high convective Mach numbers, shown in Figs. 12-16. In these cases, 
the increased compressibility results in steep pressure gradients and in the creation 
of “eddy shocklets.” These shocklets are initiated at the shear zone of the layer, and 
extend to the outer region of the flow near the boundaries. The layer is dominated by 
regions of locally subsonic and supersonic flow. To adjust to the pressure differences 
between these regions, a shocklet, albeit a weak one, is created. Further examination 
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shows that the regions of high compressibility tend to push the shocklets backwards 
in the opposite direction of the motion of the large scale structures. Therefore, the 
shocklets rotate in a direction opposite to the rotation of the large scale vortices. 
Contours of the instantaneous Mach number show that the shocklets are related to 
the formation of large scale vortical structures, illustrated in Figs. 17 through 20. The 
form of these structures, as well as localized regions of subsonic and supersonic flow 
are clearly seen. Examination of the instantaneous Mach number contours indicates 
that the strength of the shocks increase with the increase of the convective Mach 
number. 

The appearance of shocklets in flows of M c > 0.7 has been previously reported by 
Lele [9] and by Sandham and Reynolds [10] for two dimensional simulations. This 
phenomenon has not been observed in experimental studies, nor has it been reported 
in existing three dimensional mixing layer simulations. Shocklets have also been noted 
in simulations of homogeneous turbulence [29]. It has been suggested [30] that the 
appearance and strength of shocklets within the flow might be predicted by examining 
the root mean square of the Mach number and the normalized root mean square of 
the density. For example, the ratio of the two parameters may be defined as 



For Moo 0(1), if a > 1 the flow would be quasi-incompressible, i. e. no shocklets 
would appear. However, shocklets would appear for a < 1. Profiles of a are shown in 
Fig. 21 for M c = 0.2 and in Fig. 22 for M c = 0.8. The figures show that a > 1 for all 
y* when M c = 0.2, and a < 1 at some locations for M c = 0.8. Shocklets appear for 
M c = 0.8, but not for M c = 0.2. This indicates some correlation between the order of 
ratio of the root mean squares of the Mach number and density and the appearance 
of shocklets in the flow. However, more work is clearly needed. 


20 




3.2 Effects of Reaction Exothermicity 


The effects of an exothermic chemical reaction can be depicted by repeating the 
procedure above for varying values of the heat release parameter, C e , while keeping 
other parameters constant. Results are presented for constant rate kinetics with heat 
release values of C e = 0, 1.5, and 6; and for Arrhenius type k^eti'-.n w’ +Vi Ce = 1 5 **id 
Ze = 10. In these simulations, a grid of 127 x 127 points was used, the Damkohler 
number was set equal to 10, and M c was set equal to 0.2. No harmonic forcing was 
added. 

The results of these simulations are first presented in the form of plots of the 
vorticity thickness versus normalized time for all four cases (Fig. 23). This figure 
shows that the rate of growth is highest for the non-heat releasing case ( C e = 0), 
and as the heat release parameter is increased, the coupling between the reaction and 
the hydrodynamics causes the layer to grow at a lower rate. The relatively smooth 
regions of the vorticity thickness growth may be attributed to diffusion thickening 
and a jump in magnitude of this thickness represents vortex roll-up. For the C e = 0 
case, the layer responds to perturbations fairly quickly, and vortical structures are 
formed at t" % 3. An increase in the magnitude of the heat release results in a delay 
of vortex roll-up, and the jump in vorticity thickness does not occur until t * ss 7. 
Further increase in the magnitude of the heat release results in additional delays, as 
can be seen from the case of C e = 6. This is also observed for the Arrhenius model 
with C e = 1.5. In these two cases, the effects of exothermicity is most pronounced; 
vorticity roll-up does not occur at all, and the only growth in the thickness of the 
mixing layer is due to molecular diffusion. 

The vorticity contour plots demonstrate this point further. These are shown in 
Figs. 24 through 33 for each value of the heat release parameter at various times. As 
mentioned above, the non-heat releasing layer goes through roll-up and pairing fairly 
quickly (Figs. 24-26). The resulting large scale structure then rotates clockwise as 
time progresses. After the collapse of the vortex, no additional roll-up occurs, causing 
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a fluctuation in the magnitude of the vorticity thicknesses as shown in Fig. 23. The 
vorticity contours of the C e = 1.5 and constant chemical kinetics case for several 
times are presented in Figs. 27-29 . At time t* = 6, when the non-heat releasing layer 
has already rolled up, this mixing layer is only showing the initial stages of instability. 
When the large scale structures are finally formed, the vorticity thickness of the layer 
has grown via diffusion to twice that of the non-heat releasing case rigui, before die 
roll-up. As in the previous case, there are two distinct structures in the flow which 
combine into one. This is not apparent from the vorticity thickness profile, since 
the initial pairing is masked by the size of the shear layer. Further increase in the 
heat release prevent the mixing layer from responding to background perturbations 
at all. The contour plots of the vorticity in both cases, shown for t* = 8 in Figs. 31 
and 33, are composed of parallel lines, which indicate the lack of formation of any 
vortical structures. At larger times, the layer becomes too “thick” to respond to the 
background perturbations, and proceeding in time does not produce any substantial 
enhancement in mixing except that facilitated by diffusion. The mixing rate in the 
reacting layer with the Arrhenius model is severely retarded compared to that of the 
constant kinetics case. With the application of the Arrhenius reaction model, the 
rate of increase of temperature is substantially more than that of the constant rate 
kinetics model, even though the heat release parameter is kept fixed. This increase 
in the local and the global magnitudes of the temperature stabilizes the flow, and 
for the case of Ze = 10, the instability modes in the layer do not seem to grow fast 
enough to form large scale coherent vortices. 

The influence of the heat release on the structure of the flame is demonstrated by 
examining the product thickness of the layer, defined by the normalized total product 
concentration of the layer as a function of time (Fig. 34). An examination of this figure 
shows the influence of heat release on all stages of the development of the layer. At 
initial times, the effect of heat release is an enhanced product formation, while the 
reverse applies at the intermediate and final stages. Initially, the effect of heat release 
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is to expand the fluid at the cores of the layer. A large mixing zone is formed, which 
results in an enhanced reaction and an increased product formation. This explains 
the increased thickness at initial times. However, as the heat release increases and the 
layer thickens, the growth of the instability modes become subdued, postponing the 
formation of large scale vortices. After initial times, the non-heat release (C e = 0) and 
C e = 1.5 cases predict a sharp increase in the product thickness. This is caused ! Uj 
the dynamics of the large scale structure formation. The mechanism of roll-up causes 
the reactants to be entrained into the layer from their respective free streams. This 
produces an increase in the extent of mixing, reaction and product formations. The 
lack of roll-up in the C e = 6 and Arrhenius cases means the only mixing mechanism 
is due to molecular diffusion. Since in these simulations the diffusion mechanism is 

not as efficient as the roll-up of vortices in mixing reactants, the product thickness is 
correspondingly lower. 

Further influences of heat release become apparent by examining the effect on 
statistical quantities. In Figs. 35 through 37, the normalized profiles of the mean 
streamwise velocity component are presented at three different times. At early times, 
the gradients of the velocity for the low heat release cases are steeper than that of 
the higher heat release cases. This is caused by the local expansion of fluid inside the 
core of the layer. At later times the high heat release cases have a higher profile gra- 
dient, and thus, less mixing. This has a substantial influence on the two-dimensional 
turbulence transport, as indicated by the cross-stream variations of the streamwise 
velocity mean square, presented in Figs. 38 through 40. These figures show how the 
amplitude of the fluctuations decrease when the heat release is increased. Similar 
trends are observed in the profiles of the Reynolds stress shown in Figs. 41 through 
43 for three different times. The heat release clearly has a stabilizing effect on the 
flow, demonstrated by the reduction of the turbulent fluctuations. 

The results of these simulations are consistent with the linear stability analysis of 
a stratified mixing layer, in that the density reduction in the middle region of the layer 
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has a stabilizing effect [31]. However, this analysis was presented for a non- reacting 
layer. The effects of chemical reaction on linear stability is currently being studied 
by Hu and Givi [32]. 

3.3 Flame Extinction 

The nonequilibrium effects leading to local flame extinction in the reacting mixing 
layer was the final topic of this investigation. Simulations were performed using an 
Arrhenius chemical kinetics model, and repeated under otherwise identical conditions 
using constant rate kinetics for comparison. The Zeldovich number was chosen as 20, 
the Damkohler number was set to 10, the convective Mach number was set equal 
to 0.2, and the heat release parameter was set at 1.5. A grid of 127 x 127 points 
was used. The temperature field was initialized with a Gaussian distribution in the 
transverse direction and a constant distribution in the streamwise direction. The 
maximum temperature occurred along the centerline of the layer and had a value of 
five times the free stream temperature. 

The initial temperature distribution has a stabilizing effect on the flow. In order 
to trigger the formation of large scale structures, explicit harmonic forcing was added. 
The effects of three initial perturbation mechanisms were examined. These were: case 
I, the most unstable mode was added to the initial velocity profile, case II, the most 
unstable mode plus its first subharmonic were added, and lastly a control case with 
no forcing. The amplitude of the forcing, t\ for the most unstable mode, and e 2 for 
the first subharmonic of the most unstable mode, were set to 0.5% of the mean flow. 
For the purpose of evaluating the effects of harmonic forcing, the Damkohler number 
was set to zero. 

The influence forcing has on the development of the vorticity thickness is shown in 
Fig. 44. As expected, the non-forced layer remains stable and only grows via molecular 
diffusion. The forced cases appear to have identical growth until time t* ^ 1. At this 
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time, the growth of the layers follow different paths, with the case I layer growing 
to about twice the thickness of the case II and non-forced layers. This is explained 
by presenting the vorticity contours for three points in time; t* = 1, right before 
the layers diverge, t* — 1.5, soon after the split, and t* = 2.1, at a time later than 
the divergence. For the unforced case, the vorticity is shown as parallel lines; there 
is no significant instability within flow (Figs. 45 and 46). As time progresses, the 
layer thickens via molecular diffusion until there is no possibility for the background 
perturbations to grow. At t* = 1 , the forced layers, shown in Figs. 47 and 50, have 
both rolled up into a pair of vortices. It can be seen in Fig. 48 that at t* = 1.5, the 
case I vortices have elongated slightly and have taken on an elliptical shape. In case 
II, the addition of the first subharmonic initiates a second roll-up (Fig. 51) resulting 
in a pairing of the vortices. At time t* = 2.1 the case I layer, shown in Fig. 49, has 
taken an even more exaggerated elliptical shape than that at the previous time. The 
case II flow has paired into a single coherent structure shown in Fig. 52. In the center 
of the structure, the cores of the previous vortices may still be seen. 

After the examinations of the effects of harmonic forcing, simulations were then 
performed for the Arrhenius and constant rate kinetics models. Since the most pro- 
nounced mixing is desired, perturbations associated with the most unstable mode plus 
the first subharmonic were used in these simulations. Contour plots of the instan- 
taneous reaction rate and the product concentrations may be examined to provide 
insight into the progress and the structure of the reaction. These contours are shown 
in Figs. 53 and 54 for time t* = 1 , when the vortices have just been formed. Note 
the high concentrations of product in the core of the vortices where the most mixing 
occurs. The reaction rate is highest along the mixing surface of the layer, that is, in 
the center of the intervortex braids. The contours at t* = 1.5 (presented in Figs. 55 
and 56) portray the behavior at the initial stages of pairing. The reaction rate has 
begun to decrease in the braids of the emerged vortex. At time t* = 2.25, shown in 
Figs. 57 and 58, the layer has completed the pairing process and strong gradients are 
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observed in the vortex braids. Since the braids are “stretched ” as the vortex rotates, 
they are the area of the highest strain. The flame at the onset of extinction is shown 
in Figs. 59 and 60. Note that the reaction rate has gone to zero at the braids. The 
product concentration contour shows that no product exists in these extinguished re- 
gions. This demonstrates that the flame did not quench due to depletion of reactants. 
A* <,hia point in time, the flame is not continuous, forming what will be called for lack 
of a better term, a “flame eddy.” 

This extinction phenomenon has been explained by Peters [33] as follows: At the 
regions of high strain, the reactants are supplied at a faster rate than they can be 
consumed by the flame. Thus the local temperature in that area drops and the flame 
becomes very rich with the reactants. As a result, the flame is quenched in that area. 
If a fast chemistry model or an equilibrium chemistry model is used, the extinction 
mechanism can not be captured. Investigation of such phenomena requires finite-rate 
chemistry simulation in the form presented here. 
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Chapter 4 


Summary and Conclusions 


This work deals with direct numerical simulation of a compressible, temporally de- 
veloping, reacting plane mixing layer. Several simplifying approximations were made 
so that the effects of the variation of isolated parameters could be studied in de- 
tail. In particular, the chemical reaction was assumed to be of the type A + B — > 
Products -f- Heat , and thermodynamic properties were assumed constant and iden- 
tical for all the species. Two types of kinetics models were used, one with constant 
rate kinetics, and another with an Arrhenius model. In the constant rate case, the 
Damkohler number was the parameter that described the reaction, while both the 
Damkohler number and the Zeldovich number were used to describe the Arrhenius 
reactions. A simple linear gradient model was employed to model all the diffusion 
processes. Integration of the governing transport equations was performed by higher 
order finite difference methods. The two methods used here were the Gottlieb- Turkel 
two-four dissipative scheme and Carpenter’s dissipative compact parameter scheme 
(DCPS). The results obtained were not affected significantly by the choice of numer- 
ical algorithm. 

Studies of various flow phenomena were performed by varying one representa- 
tive nondimensionalized parameter while keeping all other parameters constant. The 
convective Mach number, M c , was used to describe compressibility, the heat release 
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parameter, (7 C , represented the exothermicity of the reaction, and the Damkohler 
number and the Zeldovich number quantified the extent of the reaction. 

The simulations concerning the effect of compressibility on the mixing layer showed 
a direct correlation between increased compressibility, increased stability, and reduced 
turbulence. When the convective Mach number was increased, the rate of growth of 
the mixing layer, widen <vas measured by the growth of the vorticity thickness, was 
markedly reduced. Degradation of the development of the streamwise fluctuating 
velocity and the Reynolds stress profiles was also noted. For high compressibility 
cases, eddy shocklets were observed within the flow. This has been reported in some 
previous two dimensional simulations, but has not been observed in experiments or 
in three dimensional simulations. An expansion of this problem into three dimensions 
is suggested for future work in this area. 

Increased exothermicity was observed to slow the growth of large scale structures. 
At the initial stages of development, high heat release increased the amount of prod- 
uct formed via volumetric expansion of the core of the layer. However, the heat 
release caused the layer to be less responsive to perturbations, reducing the growth 
of the layer at later stages. The overall effect of increased heat release, therefore, 
was to stabilize the flow and to decrease the extent of reaction. This was shown by 
examining the contour plots of selected quantities as well as statistical variations of 
those quantities. 

The selection of chemical kinetics model was shown to have significant effect on 
the development of the flow. The introduction of an Arrhenius chemistry model had 
a stabilizing effect, thus degrading the progress of the reaction. When the layer was 
harmonically forced, the structure of the flow was found to be controllable by varying 
the type of perturbation. When only the most unstable frequency of the layer was 
added to the mean flow, the layer went through a single roll-up process. The inclusion 
of the first subharmonic of that frequency caused the layer to go through a second 
roll-up, in the form of pairing of the neighboring vortices. Non-equilibrium effects of 
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the kinetics model was studied for the last case. It was found that at high Zeldovich 
numbers, the flame would be quenched at regions with large local values of the strain 
rates. Suggestions for future work in this aspect of research include extending the 
simulations to three dimensions and adding more realistic chemistry models. 
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Figures 



Figure 1: Schematic diagram of a temporally evolving mixing layer. 
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Figure 2: A representative contour plot with the computational grid superimposed. 
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Figure 3: Normalized vorticity thickness versus normalized time for different values 
of the convective Mach number. 
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Figure 4: Profiles of normalized mean streamwise velocity for different values of the 
convective Mach number at time t * = 6. 
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Figure 5: Profiles of normalized mean streamwise velocity for different values of the 
convective Mach number at time t m = 8. 
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Figure 6: Profiles of normalized mean streamwise velocity for different values of the 
convective Mach number at time t * = 10. 



Figure 7: Profiles of normalized mean squared streamwise velocity for different values 
of the convective Mach number at time t * = 6. 
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Figure 8: Profiles of normalized mean squared streamwise velocity for different values 
of the convective Mach number at time f* = 8. 



Figure 9: Profiles of normalized Reynolds stress for different values of the convective 
Mach number at time t* = 6. 
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Figure 21: Profile of ratio of the rms Mach number to the normalized rms density for 
M c = 0.2 
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Figure 22: Profile of ratio of the rms Mach number to the normalized rms density for 
M c = 0.8 
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Figure 23: Normalized vorticity thickness versus normalized time for different values 
of the heat release parameter. 
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Figure 24: Plot of vorticity contours at time <* = 6, C t = 0 



Figure 25: Plot of vorticity contours at time t* = 8, C e = 0 
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Figure 28: Plot of vorticity contours at time t" = 8, C e = 1.5 



Figure 29: Plot of vorticity contours at time t’ = 10, C e = 1.5 
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Figure 30: Plot of vorticity contours at time t m = 6, C e = 6 
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Figure 31: Plot of vorticity contours at time t’ = 8, C e = 6 
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Figure 33: Plot of vorticity contours at time t m ~ 8, C e = 1.5, and Arrhenius kinetics 
model. 
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Figure 34: Normalized product thickness versus normalized time for various values 
of the heat release parameter. 



Figure 35: Profiles of normalized mean streamwise velocity for different values of the 
heat release parameter at time t m = 3. 
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Figure 36: Profiles of normalized mean streamwise velocity for different values of the 
heat release parameter at time t m = 6. 



Figure 37: Profiles of normalized mean streamwise velocity for different values of the 
heat release parameter at time t m = 8. 






Figure 38: Profiles of normalized mean squared streamwise velocity for different values 
of the heat release parameter at time t” = 6. 
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Figure 39: Profiles of normalized mean squared streamwise velocity for different values 
of the heat release parameter at time t* = 8. 
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Figure 40: Profiles of normalized mean squared streamwise velocity for different values 
of the heat release parameter at time t' = 10. 
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Figure 41: Profiles of normalized Reynolds stress for different values of the heat 
release parameter at time t* = 6. 
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Figure 42: Profiles of normalized Reynolds stress for different values of the heat 
release parameter at time t m = 8. 
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Figure 43: Profiles of normalized Reynolds stress for different values of the heat 
release parameter at time t x = 10. 
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